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Abstract 

Wc develop a technique to compute the high-frequency asymptotics of spin cor- 
relators in weakly interacting disordered spin systems. We show that the dynamical 
spin correlator decreases exponentially at high frequencies, (SS)^ ^ exp(— r*a;) and 
compute the characteristic time r* of this dependence. In a typical random configu- 
ration, some fraction of spins form strongly coupled pairs, which behave as two-level 
systems. Their switching dynamics is driven by the high-frequency noise from the sur- 
rounding spins, resulting in low-frequency 1// noise in the magnetic susceptibility and 
other physical quantities. We discuss application of these results to the problem of 
susceptibility and flux noise in superconducting circuits at mK temperatures. 



Introduction 

In many physical systems, the relevant degrees of freedom are discrete degenerate quantum 
variables such as spins (electron or nuclear) that interact weakly with each other. At suf- 
ficiently high temperatures, the spins are completely disordered and characterized by the 
trivial thermodynamic correlators. However, their dynamical properties remain interesting 
and important for various applications. The spins generally fluctuate at frequencies of the 
order of the typical interaction Jtyp. For some purposes, though, one needs to know the 
asymptotics of the spin correlation function at much higher frequencies. In this regime, the 
fluctuations and dynamical response are due to the simultaneous excitation of many degrees 
of freedom and inherent nonlinearity of the system. The goal of this paper is to develop a 
formalism for the computation of the high-frequency spin correlator and to apply the result 
to the problem of magnetic (flux and inductance) noise at mK temperatures produced by 
paramagnetic spins located at the surface of superconductor. 

Speciflcally, we consider a set of spins described by the Heisenberg Hamiltonian with 
random couphngs at high temperatures, T » Jtyp. We show that the spin correlator decays 



exponentially with frequency at a; ^ Jtyp, i.e. (SS)^ ~ exp(— cw/ Jtyp), and compute the 
coefficient c in this formula. 

The smallness of the high frequency fluctuations has an important consequence when the 
interaction, although of the order of Jtyp for the typical pairs of spins, is much larger for a 
small but significant number of pairs. This happens, for instance, if the spins are positioned 
randomly while the interaction between them falls off fast with distance. Each strongly 
coupled pair has two states, the singlet and the triplet. Transitions between those states 
require the exchange of a large energy quantum, AE = J ^ Jtyp between the pair and the 
surrounding spins. The transition rate is proportional to the dynamical spin correlator at 
frequency AE, which is exponentially small. As a result, the switching events are rare and 
produce noise with, as we show below, 1/f spectrum at very low frequencies. At intermediate 
frequencies, the 1/f noise is due to spin diffusion [1]. 

Our arguments rely on two important assumptions. First, the condition on the tem- 
perature is stricter than in the high frequency case: we need that T ^ Jtyp lii( Jtyp/i^) if 
the noise is measured at frequency u. The right-hand side of this inequality represents the 
transition energy AE of the relevant two-level systems; if it is too large, such systems are 
frozen. We also assume that there is no efficient energy exchange mechanism between the 
two-level systems and an external thermal bath, e.g. electrons or phonons. 

One application of our theory is the property of paramagnetic spins on the surface of 
superconducting aluminum. These spins, presumably located in the thin AI2O3 film or just 
at the metal-insulator interface, are believed to be responsible for the low frequency fiux 
noise that limits the sensitivity of dc SQUIDs [21 [3] and causes decoherence in superconduct- 
ing fiux qubits |1H6]. It is very likely that the coupling between the spins is due to RKKY 
interaction [T], so their dynamics is controlled by the Heisenberg Hamiltonian. A number 
of experiments j3l [71 [8] shows that the 1/f noise at very low frequencies is temperature- 
independent fot T ~ 50-500 mK. We take these data as evidence that the inter-spin inter- 
action Jtyp is so small that the spins remain in the high temperature regime in the whole 
experimental range. Finally, in the superconductor all electronic excitations are gapped, 
and direct interaction between the spins and phonons at low temperature is negligible. All 
these facts combined together indicate that the spins in the oxide layer on the surface of 
aluminium or similar superconductors are described by the model studied in this paper. 

A lower bound on Jtyp can be inferred from the recent measurements of the fiux qubit 
relaxation rate [9j. If we assume that the relaxation is mainly due to the interaction with the 
surface spin, then Freiax oc Im(S'5')^=A£;; where AE is the energy difference between the two 
qubit states. The experimental data indicate the persistence of magnetic relaxation at very 
high frequencies, at least up to 1 GHz. We conclude that the typical interaction between the 
surface spins is Jtyp > 1 GHz ~ 50 mK. This is consistent with the observed spin freezing at 
~ 50 mK in some samples [Tj as well as a theoretical estimate of the RKKY interaction, see 
Section IH 

The theory developed in this paper has many other applications, besides the problem of 
magnetic noise in superconducting circuits. It can be applied to many problems in which the 
relaxation of high energy modes is possible only by the simultaneous creation of a very large 
number of low energy excitations. Some examples are the relaxation of inverted nuclear spin 
polarization in high magnetic fields and a similar problem of magnetization relaxation of cold 
atoms following a rapid quench |10| or 1// magnetic noise observed in a low temperature 
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calorimeter with superconducting particle absorbers [TT]. 

The paper is organized as follows: we first study high frequency spin fluctuations at 
infinite temperature; then we show how the high frequency noise generates low frequency 
noise in a model with a broad distribution of coupling parameters. We finally discuss the 
application of this model to the problem of flux noise in superconducting qubits. 



1 High frequency noise 

We consider a set of spins of magnitude S described by the Heisenberg Hamiltonian, 

H = —'^^JjkSjSk (1) 

j<k 

with random couplings Jj^. We are primarily interested in the S = 1/2 case because it is most 
relevant for the physical applications. Let us fix j and consider the normalized correlation 
function for the spin Sj (with components Sj, a = 1,2, 3) in the time and frequency domains: 

Fit) = j^^^ {S]{t)S]{0)) , F{uj) = j F{ty-Ut. (2) 

Note that F{ijj) is real and positive. This follows from its physical interpretation as the 
quantum transition probability under a small perturbation V{t) = —h°'e~^'^^Sj + h.c. (More 

abstractly, Jj|^(X^(t)X) (it > if the averaging is performed over any quantum state that 
commutes with the Hamiltonian.) 

In general, the asymptotics of the Fourier transform at large u are related to analytic 
properties of the original function. However, since F{u}) > 0, it is sufficient to consider F{t) 
on the imaginary axis: F^ir) = F{u)e'^'^^. If the integral converges for < r < r*, then 

F{u) should decay sufficiently fast, and the real-time correlator has an analytic continuation 
to the strip < Imt < r^,. If, in addition, the integral diverges for r > r*, then F[u}) ~ e~'^*'^ 
for u — )■ +00. (Strictly speaking, this expression is valid if we smooth out some unimportant 
features, such as peaks at discrete frequencies.) 

At infinite temperature, the operator average in ^ is taken over the maximally mixed 
state, i.e. (X) = TrX/ Tr U. The cyclic property of the trace implies that both variants of 
the correlation function are even and real. Therefore, the Taylor expansion of the real-time 
spin correlator contains only even terms: 

m = (e-*^;(0)e--^;(0)> = f: ^^H^t-, (3) 

■ir / \ /•+~ ^, ^ 



We are interested in the asymptotic behavior of the moments for n — t- oo. Let us assume 
that the power series ^ has a finite convergence radius r*. Then F{t) has singularities at 
t = ±ir^,, and the frequency representation of the correlation function decays as e"'^*''^' for 
large u. Our main goal is to show that r^, is nonzero and finite, and to compute its value. 
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Physically, the spin correlation function does not completely vanish at high frequencies 
only due to nonlinearities that allow for combining the precession frequencies of a large 
number of spins. To understand it qualitatively, consider an infinite system in which each 
spin is coupled to Z neighbors with interaction strength Jtyp. If an energy quantum lo ^ Jtyp 
is exchanged between a single spin and some external object, it has to be distributed between 
other spins and divided into portions of the order of Jtyp. The most efficient processes of 
this kind are cascades that correspond to tree-like diagrams and do not involve any back 
reaction. Summing over such diagrams is equivalent to the mean-field approximation. It 
is generally applicable for Z ^ 1, but will also produce qualitatively correct results if the 
energy has many paths to escape, for example, on a tree with Z > 2. 

Thus, we assume that each spin j precesses in the local exchange field hj(t) = Ylk JjkSk{t) 
produced by the spins surrounding it. The fiuctuating field h{t) = hj{t) is characterized by 
the correlation function 

{h'^{t)h\0)) = g'F{t)6^\ (5) 



where a, b are spatial indices (a, b = 1,2, 3) and 



^^^i^'^k J'jk- Using the mean-field 



approach, we assume that h{t) is Gaussian. The relevant energy scale in this model (corre- 
sponding to Jtyp in the qualitative discussion above) is given by g. In this approximation, 
the time evolution of each spin simplifies to dS /dt = S x h{t), or, in a matrix form. 



dsydt = [h{t)us'' 



(6) 



Here h is the matrix corresponding to the cross-product with h, and [h]ca are its elements: 

/ -h^ \ 

h=[ -h^ ], [hU = ecabh'. (7) 

V -h' / 

We solve equation (6) in a symbolic form, S^{t) = Texp h{t) dt ^"(O) and average over 

I I L J ca 

the thermal state of the spin: {S^S^) = ^^^6ac. Thus, the spin correlator is expressed as 
an average over the fiuctuating fields: 



F{t) = ^ {Tr (^Texp J h{t)dt 



(8) 



Using equations ([s]) and (jsj), we recursively compute the first few moments D2n- 



Do = l, 

D2 = 2g'Do, 

D, = 2g^D2 + lOg^Dl 

= 2g^D^ + ASg^D^Do + rOg^Dl 



(9) 



It is unfortunately difficult to pursue this computation much further and find the asymptotic 
behavior of the moments. Nevertheless, these recursive formulas give an exact expansion of 
the spin correlation function in lower orders: 



F{t) = l-igtr + ^igtr 



97 
360 



{gtf + Oigtr. 



(10) 
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Let us describe a general method for the calculation of the moments. It can be used to 
derive equations ^ and continue to higher orders, though we could not obtain a closed-form 
expression for D2n- We first use the equation of motion (|6| to express the 2?T,-th derivative 
in Eq. Q: 



2n 

E 

k=l 



2n-l 
k-1 



where m'^^ = dl^h/dt^. By applying Eq. (11) repeatedly, we obtain an expansion 



dt 



2n 



m>2 ki,...,km>l 
fclH hkm=2n 



with the coefficients 



(2n- 1)! 



{2n-h)---{2n-h km-i) (h - 1)\ ■ ■ ■ {k^ - 1)1' 



(11) 



(12) 



(13) 



Recall that the spin operators in ( 12 ) pertain to site j, and the fields h are linear combinations 
of adjacent spins. The field derivatives h^''^ can be expressed in terms of further neighbors. 
So far all the calculations have been exact. Now we employ the mean- field approximation 
and average over the thermal state of the spin, and then over the Gaussian fields. Thus, 



D 



2n 



E E (Tr(/.('=-i)/.(^-i) ■ ■ • /.(^-i)) 



m>2 A;i,...,fcm>l 
fciH \-km=2n 



(14) 



The Gaussian average in this formula can be evaluated using Wick's theorem and the explicit 
form of the second-order correlator: 



(15) 



2 Estimates for r* 



The difficulty in carrying out the above calculation lies in the matrix nature of the fields h. 
We now introduce two different approximations, which are expected to shift the answer in 
opposite directions by overestimating or underestimating the moments D2n, thus giving a 
lower and an upper bound for r^,. Remarkably, we find that the two bounds are quite close 
in value. 

In our first approximation, we ignore the noncommutativity of the matrices h^''\ or 
equivalently, h{t) for different values of t. More exactly, we symmetrize over all permutations 



of those matrices. Let us briefly explain how to use this recipe for equation (14), although it 



is easier to start almost from scratch as we will do later. By symmetrizing the combinatorial 
coefficients (13), we obtain the expression , ^.^"^^ , . Then we consider an individual Wick 
pairing for the product /^^^i^^) ■ ■ ■ h^^rn-i) g^j-^^ relate it to the corresponding paring for the 

k~l J 

m-th power of h. The ratio is given by the scalar coefficient D^^i in Eq. (15|); there 



is one such factor for each pair. Summing over all Wick pairings, we get (Tr /;,"^times a 
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certain number. It is easy to show that (Tr/i^^') = {-g'^Y^^^ for r > 0. Eventually, 
Eq. (14) becomes: 

An _ 2 ^ (2r + 1) ^ ^Dg^ 

(2n)! 3 ^ r! ^ ^ gi! g,! ' ^ ^ 

giH l-gr=2n 

The solution to this recurrence yields the following expansion of F{t) in lower orders: 

7 qq 

F(t) = 1 - {gty + -{gtr - —{gtr + ... (17) 



Notice the subtle difference from the exact solution (10). We believe that also in higher 
orders, the noncommutativity of the field matrices results in smaller values of than the 
ones obtained by the present method. 

The same approximation can be applied directly to Eq. ([s]) if we replace the time ordered 
exponential with an ordinary exponential: 

F{t) = ^ (Tr(exp^(t))) = ^ (l + 2cos|#(t)|) , H{t) = j\{t) dt, (18) 

where H is the matrix associated with vector H. The vector components H"-{t) (a = 1, 2, 3) 
are Gaussian fields with the following second moment: 

{H''{t)H\t)) = 2g^G{t)5''\ (19) 

G{t) = [ {t-ti)F{ti)dti. (20) 
Jo 

Higher moments {\H\^^) (for a given t) are calculated directly by integrating over H . Using 
the notation x = \H\/{2g'^G{t)Y/^ , we find that 

,2n. ^ jTx'-e-^'/'xMx ^ r(n + 3/2) ^ (2n + l)! 

^ ' j^e-V^xMx r(3/2) 2"n! ' ^ ' 



Now we calculate the average in Eq. (18) by expanding cos\H{t)\ in powers of \H\^. The 
result is: 

F(t) = ^ (l + 2 (1 - 2/G(t)) e-^'«W) . (22) 

This equation needs to be solved in conjunction with Eq. (20). 

To compute F{it), we formally replace t with r and —g^ with +(?^. The numerical 
solution exhibits singularities at r = ±r* for gr^, ^ 1.74. Since the moments D2n have been 
overestimated, this gives a lower bound for the true value of r*. 

We now make a different approximation to matrix products, exaggerating the noncom- 
mutativity. Specifically, we replace the physical = 1/2 spin for the SU{2) group with an 
SU{N) spin and consider the limit — t- oo. Let us first introduce the generators of U{N): 
these are N x N matrices cr"^ («,/? = 1, . . . ,N) with the elements [o""'']^,^ = 5a^5py. The 
SU (N) spin is defined by the traceless matrices 

= _L^»I3^_ (23) 
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b) 



^ ] yfi J |iv ~ 




d) 



Figure 1: (a) A generic term in the diagrammatic expansion of (TT(uj^{t) a"'^ UN{t) a'^"')) 
(for an arbitrary A^) using equations (25) and (26). The arrows indicate f/Ar(t) (the propa- 
gation from to t) and Ulf{t) (moving backward in time); the shaded ribbons correspond to 
field correlators, (b) The exact definition of lines and ribbons. Note that that they do not 
depend on time: the factors F{tk — tj) are coefficients in front of the diagram, (c) A planar 
diagram: one of the leading terms in the A^ — > oo limit, (d) A diagram for the function f{t) 



defined in (28). 



The generalized Heisenberg interaction is given by the operator ^cr^'^af" with the matrix 

elements ;^ [cr"^] [a'^"'] = -^(^S^j,6s^ — jf^'ysS^u^. Note that if we let a"^ and 5"'^" act on 

the same spin (by contracting over a pair of indices, 5 = /i), we get 1 — ^ times the identity 
matrix. For N = 2, this matches the value of S°'S"' = S{S + 1). 
Thus, the full Hamiltonian has the form 



j<k 



(24) 



Each spin j experiences the effective field hj{t) = Ylk Jjk^^t ■ mean-field approx- 

imation, we have: 

(r^(t)/i'^'°'(0)) 



N 



FN{t) S' 



- —5" 
N 



(25) 



where = Jff^, and the spin correlator is normalized so that Fn(0) = 1. Note that 

the last term in (25) has no effect when we consider the action of the field on the spin, 
because ^"/^a^" = 0. 
The calculation of 

FN{t) = j;^^ (TT{uj,{t)d''^UN{t)a^'')), f/^(t) =Texp l'\h^%h)a^^dh (26) 



can be performed using the standard diagrammatic technique for matrix models, see Fig.[T]a,b. 
The ribbons in the diagrams may actually be changed to double lines. To justify this step. 
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let us express cx"^ in terms of cr"^ using g). The result is: (Tr a^^^ f/^(t) ) 



(Tr(Ul^[t) a"'^ U]sr{t) a^") ) ~ which shows that the ribbon between and t can be changed 
as stated (up to the —1 term). We can likewise replace the a^'^ that enter the expressions for 
UN(t), Ulfit). The contributions of the unit operator in (23) to UN(t) and to Ul^it) cancel 



each other, therefore the remaining ribbons are also equivalent to double lines. 

A great further simplification occurs in the limit of large |15) . The leading contribution 
to F]\f{t) comes from diagrams that maximize the number of loops for a given expansion order 
in g'^/N. Such diagrams are planar, and the time evolution factors UN{t), Uj^{t) separate; 
an example is shown in Fig. [l]c. Thus, in the limit N 00 equations (25) and (26) become: 



Fooit) = fit), f{t)l = (Texpj\h^^{T)a^^dT^ . 



(27) 
(28) 



The "spin propagator" f{t) can be computed by summing up planar diagrams. Specifically, 
f{t) is represented by a single line going from to t and dressed with double lines (the field 
correlators) attached on one side; these lines do not intersect, see Fig. [T]d. This diagrammatic 
expression leads to the Dyson equation: 



fit) = 1-9' 



dt2dtj%t2 - ti) /(ti 



(29) 



t>t2>ti>0 



Curiously, the Taylor expansion for -Foo(t) = f (t) coincides with the result (10) for F2{t) 
up to the sixth order. The numerical solution of Eq. (29) exhibits singularities at t = ibir* 
for gr^: ^ 1.78. Thus, the exact value of r^, for the original problem should be between 1.74/ g 
and 1.78/g. Recall that this number is the exponent in the expression for the high frequency 
noise spectrum: (/i^) oc e~'^*'^'. 



3 Low frequency noise 

The high frequency fluctuations of local effective fields drive the switching dynamics of 
strongly coupled spin pairs, and thus determine the noise at low frequencies (provided the 
number of such pairs is significant). We will find the spectrum of the low frequency noise in 
a specific model, which has been previously discussed in [Ij. 

Let us consider randomly distributed spins with the two-dimensional density p2D on the 
surface of a superconductor. The coupling of each spin to the conduction electrons is de- 
scribed by the Kondo term in the Hamiltonian: J'j S^s'^I'Tj), where s"(r*) = | [(7°]^,^ ^)t(r)^i^(^) 
The relative coupling strength can be characterized by the dimensionless parameter 

X = Ju={ln{Ep/TK)y\ (30) 

where u is the electron density of states per spin projection and Tk is the Kondo temperature. 
The RKKY interaction between two spins can be expressed as follows: 

Jjk = - ^ 3 cos{2kFrjk) = J (rjk) cos (pjk, J{r) = - — r"^, (31) 
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where (j)jk is random, and rjk is assumed to be much smaller than the coherence length. The 
surface is usually covered by an insulator and thus can carry spins with arbitrary coupling 
to the conduction electrons, ranging from A ~ 1 down to 0. However, only the spins with 
Tk ^ Tc are active, the other being quenched by the Kondo effect. This imposes the 
constraint A < 0.1. The typical distance between the spins is rtyp = P2d'^j which defines 
the main parameter of our model: g ~ J(rtyp) = ^ • Spins coupled with this strength 
form an infinite cluster, whereas more strongly coupled pairs are rare. If we use the values 
u ^ 17eV~^nm~^ (the DOS in aluminum) and p2D ~ 0.4 nm~^ (the typical spin density 
measured in several types of films [TJ [S]), we get the estimate g < 50 mK. The upper 
bound corresponds to the spins with the largest value of A. There are some indications, both 
experimental [H] and theoretical [17], suggesting that most localized spins arise from surface 
bound state. Therefore we can expect the actual value of g to be close to the upper bound. 

The density of pairs with interaction J ^ g can be found by calculating the probability 
that a given spin has a neighbor at the corresponding distance. Such pairs constantly switch 
between the singlet and the triplet states at a characteristic rate T ~ g e~'^*'^ = ge~'^'^^^, 
where c is a numerical constant. We estimate the density of pairs in terms of F: 

^ = \p2D27rrdr ~ (J/g)-'/' diJ/g) ^ (ln(^/r))-'/'^. (32) 

Each pair generates a random telegraph noise in any quantity x that is proportional to the 
number of pairs in the triplet state; for example, x can be the magnetic susceptibility. This 
noise has a Lorentzian spectrum with width F, but, for the purpose of a crude estimate, we 
may assume that it is concentrated at frequency u = T. Therefore, the total noise spectrum 
is proportional to the previously calculated distribution of F: 

((^^)'> ^ , 5/3, / ^ - (33) 

Thus, the magnetic susceptibility noise has a roughly 1/f spectrum at low frequencies. In 
the presence of a constant magnetic field, it translates to a 1// flux noise. Some experimental 
data indicate a spontaneous breaking of the time reversal symmetry in spin systems at low 
temperatures [l_6j. Although the exact nature of the symmetry-breaking order parameter 
is unclear, it is most likely to exert local fields on the pairs of spins that undergo the slow 
fluctuations. In this case, the triplet state is split in energy and produce a net magnetization, 
observed flux noise. 
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